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BI-LEVEL INTEGRATED SYSTEM SYNTHESIS (BLISS) 


Jaroslaw Sobieszczanski-Sobieski *, Jeremy S. Agte* f 
and Robert R. Sandusky , 


Abstract 

BLISS is a method for optimization of engineering 
systems by decomposition. It separates the system 
level optimization, having a relatively small number of 
design variables, from the potentially numerous sub- 
system optimizations that may each have a large num- 
ber of local design variables. The subsystem optimiza- 
tions are autonomous and may be conducted concur- 
rently. Subsystem and system optimizations alternate, 
linked by sensitivity data, producing a design im- 
provement in each iteration. Starting from a best guess 
initial design, the method improves that design in it- 
erative cycles, each cycle comprised of two steps. In 
step one, the system level variables are frozen and the 
improvement is achieved by separate, concurrent, and 
autonomous optimizations in the local variable subdo- 
mains. In step two, further improvement is sought in 
the space of the system level variables. Optimum sen- 
sitivity data link the second step to the first. The 
method prototype was implemented using MATLAB 
and iSIGHT programming software and tested on a 
simplified, conceptual level supersonic business jet 
design, and a detailed design of an electronic device. 
Satisfactory convergence and favorable agreement with 
the benchmark results were observed. Modularity of 
the method is intended to fit the human organization 
and map well on the computing technology of concur- 
rent processing. 

0. Introduction 

Optimization of complicated engineering systems by 
decomposition is motivated by the obvious need to dis- 
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tribute the work over many people and computers to 
enable simultaneous, multidisciplinary optimization. It 
is important to partition the large undertaking into 
subtasks, each small enough to be easily understood 
and controlled by people responsible for it. This im- 
plies granting people in charge of a subtask a measure 
of authority and autonomy in the subtask execution, 
and allowing human intervention in the entire optimi- 
zation process. 

Reconciliation of the need for subtask autonomy with 
the system level challenge of “everything influences 
everything else” is difficult. Each of the leading MDO 
methods that have evolved to date (survey papers: 
Balling and Sobieszczanski-Sobieski, 1996; and Sobi- 
eszczanski-Sobieski, J., and Haftka, R. T, 1997) tries to 
address that difficulty in a different way. In the system 
optimization based on the Global Sensitivity Equations 
(GSE) (Sobieszczanski-Sobieski, J. 1990, Olds, J. 
1992, Olds, J. 1994), the partitioning applies only in 
the sensitivity analysis while optimization involves all 
the design variables simultaneously. The Concurrent 
SubSpace Optimization method provides for separate 
optimizations within the modules (Sobieszczanski- 
Sobieski, J. 1988, Renaud and Gabriele, 1991, 1993, 
and 1994; Stelmack and S. Batill, 1998) but handles all 
the design variables simultaneously in the coordination 
problem. The Collaborative Optimization method 
(Braun and Kroo, 1996; Sobieski and Kroo, 1998) also 
enables separate optimizations within the modules, 
each performed to minimize a difference between the 
state and design variables and their target values set in 
a coordination problem. This problem combines the 
system optimization with the system analysis, therefore 
its dimensionality may be quite large. 

Most of the above method implementations had to 
overcome difficulties with integration of dissimilar 
codes. This has stimulated use of Neural Nets and Re- 
sponse Surfaces as means by which subdomains in the 
design space may be explored off-line and still be rep- 
resented to the entire system. Unfortunately, effective- 
ness of this approach is limited to approximately 12 to 
20 independent variables, hence, it is best suited for the 
early design phase. Consequently, a clear need remains 
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for a method applicable in later design phases when 
the number of design variables is much larger. Meth- 
ods that build a path in design space fit that require- 
ment. Ultimately, one needs both domain -exploring 
methods and path-building methods, enhanced with 
seamless ‘gear-shifting 1 between the two. 

Motivated by the above state of affairs, BLISS attacks 
the problem by performing an explicit system behavior 
and sensitivity analysis using the GSE, autonomous 
optimizations within the subsystems performed to 
minimize each module contribution to the system ob- 
jective under the local constraints, and a coordination 
problem that engages only a relatively small number of 
the design variables that are shared by the modules. 
Solution of the coordination problem is guided by the 
derivatives of the behavior and local design variables 
with respect to the shared design variables. These de- 
rivatives may be computed in two different ways, giv- 
ing rise to two versions of BLISS. 

In either version, BLISS builds a gradient-guided path, 
alternating between the set of disjointed, modular de- 
sign subspaces and the common system-level design 
space. Each segment of that path results in an im- 
proved design so that if one starts from a feasible state, 
the feasibility in each modular design subspace is pre- 
served while the system objective is reduced. In case of 
an infeasible start, the constraint violations are reduced 
while the increase of the objective is minimized. Be- 
cause the system analysis is performed at the outset of 
each segment of the path, the process can be termi- 
nated at any time, if the budget and time limitations so 
require, with the useful information validated by the 
last system analysis. In addition to enabling complete 
human control in the subspace optimization, BLISS 
allows the engineering team to exercise judgment, at 
any point in the procedure, to intervene before com- 
mitting to the next successive pass. 

BLISS has been developed in a prototype form and has 
been successfully demonstrated on the small-scale test 
cases reported herein. 

1. Notation 

BB - black box, a module, in the mathematical model 
of a system. 

BBA(Y r ,(Z,X r )) - analysis of BB r to compute Y r for 
given Z and X r 

BBOF r - BB Objective Function computed in BB r 
BBOPT r (X r ,<t> r ,G r ) - optimization in BB r defined by 
eq.(2.1/9) 

BBOSA r (X r op{ ,Z,Y r>s ) - analysis of BB optimum for 


sensitivity to parameters 

BBSA(D(Y r ,(Z,X r ,Y r s )) - sensitivity analysis of BB r to 
compute its output derivatives w.r.t. Z, X r , and Y rs 
D(V1,V2) - total derivative dVl/dV2 
d(Vl,V2) - partial derivative 5V1/5V2; 

D(), and d() dimensionality depends on the dimen- 
sionalities of VI and V2: 

VI and V2 - are both scalars, then D and d are 
scalars 

VI vector, V2 scalar, then D and d are vectors 
VI scalar, V2 vector, then D and d are vectors 
VI vector, V2 vector, then D and R are matrices 
G 0 - vector of constraints active at the constrained 
minimum, length NG 0 

G r - vector of the constraint functions, g r t local to BB r , 
g r , t < 0 is a satisfied constraint 
G yz - constraints in a BB that have a stronger depend 
ence on Y and Z, than on X 
GSE - Global Sensitivity Equations (Sobieszczanski- 
Sobieski, 1990); GSE/OS - GSE/Optimized Sub- 
systems. 

I - identity matrix. 

L - vector of the Lagrange multipliers corresponding to 
G 0 , length NG 0 
LP - Linear Programming 
NB - the number of BBs in the system 
NLP - NonLinear Programming 
opt - subscript denotes optimized quantity 
P - vector of parameters, pi, kept constant in the proc- 
ess of finding the constrained minimum, length 
NP. 

SA((P,Z,X),Y) - system analysis; a computation that 
outputs Y for a system defined by P, Z, and X 
SOF - System Objective Function computed in one of 
the BBs 

S0PT(Z,O) - system objective optimization defined by 
eq. (2.2.3/1) 

SSA(D(Y,(Z,X)) - system sensitivity analysis to com- 
pute sensitivity of the system response Y w.r.t. Z 
and X 

TOGW - take-off gross weight 
T - superscript denotes transposition. 

X r - vector of the design variables x r j , length NX r , 
these variables are local to BB r ; X without sub- 
script - a vector of all concatenated X r , length NX 
XL, XU - lower and upper bounds on X, side- 
constraints. 

Y r - vector of behavior (state) variables output from 
BB r , these are the coupling variables; an element 
of Y r is denoted y fti ; some of y r>i are routed as in- 
puts to other BBs, and may also be routed as out- 
put to the outside; the Y r length is NY r ; Y without 
subscripts - a vector of all concatenated Y r , length 
NY 
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Y r tS - vector of variables input to BB r from BB S , these 
are the coupling variables; an element of Y r s is de- 
noted y r , s>i ; note that by this definition Y r>s is a sub- 
set of Y r , vector length NY r , s 
Z - vector of the design variables z k that are shared by 
two or more BBs, these are the system-level vari- 
ables; length NZ 

0 - subscript denotes the present state from which to 
extrapolate, or the optimal state. 

ZL, ZU - lower and upper bounds on Z, side- 
constraints 
A - increment 
AZL, AZU - move limits 

- the local objective function in BB r 
O - the system objective function equated to one, par- 
ticular yu 

2. The Algorithm 

In this section, the symbols defined in Notation are 
used in a shorthand manner without repeating their 
definitions. 



Figure 1 : System of Coupled BBs 

The algorithm is introduced using an example of a 
generic system of three BBs, as shown in Figure 1. 
Three is a number small enough for easy conceptual 
grasp and compact mathematics, yet large enough to 
unfold patterns that readily generalize to larger NB. 
Even though the system in Figure 1 is generic, it may 
be useful to bear a specific example in mind. Let it be 
an aircraft so that: 

BB 1 - performance analysis 
BB2 - aerodynamics 
BB3 - structures 

O - maximum range for given mission characteristics 
Y i 2 - includes the aerodynamic drag; Y u - includes the 
structural weight; Y 2 ,i - includes Mach number; 
Y 3( i - includes TOGW; Y 2l 3 - includes the struc- 
tural deformations that alter the aerodynamic 
shape; Y 3>2 - includes the aerodynamic loads 
g Lt - a noise abatement constraint on the mission pro- 


file; g 2<t - limit of the chordwise pressure gradient; 
g 3>t - allowable stress 

x 1( j - cruise altitude; x 2 ,j - leading edge radius; x 3 ,j - 
sheet metal thickness in the wing skin panel No. 
138 

Zj - wing sweep angle; z 2 - wing aspect ratio; z 3 - wing 
airfoil maximum depth-to-chord ratio; z 4 - loca- 
tion of the engine on the wing 

The system in Figure 1 is characterized by BB level 
design variables X, and by system-level design vari- 
ables Z. As a reference, if an all-in-one optimization 
were performed, observing the system at a single level 
and making no distinction between the treatment of X 


variables and the 
could be stated 

treatment of Z variables, the problem 

Given: 

X and Z 

(i) 

Find: 

AX and AZ 


Minimize: 

<D(X,Z,Y(X,Z» 


Satisfy: 

G(X,Z,Y(X,Z) 



Since BLISS approaches this optimization by means of 
a system decomposition, the algorithm depends on the 
availability of the derivatives of output with respect to 
input for each BB. That assumes the differentiability of 
the BB internal relationships to at least the first order. 
It is immaterial how the derivatives are computed, fi- 
nite differencing may always be used, but it is expected 
that in most cases one will utilize one of the more effi- 
cient analytical techniques (Adelman and Haftka, 
1993). 

The algorithm comprises the system analysis and sen- 
sitivity analysis, local optimizations inside of the BBs 
(that includes the BB-internal analyses), and the sys- 
tem optimization. We will not elaborate on SA beyond 
pointing out that it is highly problem-dependent, and 
likely to be iterative if there are any non-linearities in 
the BB analyses. Each pass through the BLISS proce- 
dure improves the design in two steps: first by concur- 
rent optimizations of the BBs using the design vari- 
ables X and holding Z constant; and next, by means of 
a system-level optimization that utilizes variables Z. 
We begin with the BB- level optimization. 

2.1. BB -level (discipline or subsystem) 
optimizations . 

The basis of the algorithm is the formulation of an 
objective function unique for each BB such that mini- 
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mization of that function in each BB results in the 
minimization of the system objective function. To in- 
troduce that formulation let us begin with the system 
objective function (SOF). The SOF is computed as a 
single output item in one of the BBs; without loss of 
generality we assume that it is BBj so that 

<& = yi,i (i) 

is one of the elements of Yj. 

Total derivatives of Y w.r.t. x rJ , D(Y, x rJ ), are com- 
puted according to Sobieszczanski-Sobieski, 1990, by 
solving a set of simultaneous, algebraic equations 
known as Global Sensitivity Equations, GSE, (see Ap- 
pendix, Section 1, for details) for a particular x r j 

[A] (D(Y,x r j)} = {d(Y, x r j) } (2) 

where A is a square matrix, NYxNY, composed of 
submatrices forming this pattern 


I 

A l t2 

A 1,3 

A 2,1 

I 

A 2,3 

A 3,l 

A 3 ,2 

I 


(3) 

where I stands for identity matrix, NY r xNY r , and A r s 
are matrices of the derivatives that capture sensitivity 
of the BB r output to input. For example 

A 2 ,a = -[d(Y 2 ,Y 3 )], NY 2 xNY 3 

A 3 ,2 = -[d(Y 3f Y 2 )] f NY 3 xNY 2 (4) 

One should note that eq. 2 can be efficiently solved for 
many different x r j using techniques available for linear 
equations with many right-hand sides. 

Having D(Y,x rJ ) computed from eq. 2 for all x r j, we 
can express O as a function of X by the linear part of 
the Taylor series 

^ > = yi4 = (yi,i)o+D(y u , X, ) T AX i + 

D(y 1 , i ,X 2 ) T AX 2 + D(y 1J> X 3 ) T AX 3 (5) 

where D-terms are vectors of length NX r . 

We see from eq. 5 that 

M> = D(y 1 ,,X|) T AX 1 + D(y u ,X 2 ) T AX 2 + 

D(y u ,X 3 ) T AX 3 (6) 


the three terms showing explicitly the contributions 
to AO of the local design variables from each of the 
three BBs. 

It is apparent that to minimize AO we need to charge 
each BB with the task of minimizing its own objective. 
Using BB 2 as an example, objective <j> 2 is 

♦2 = D(y u , X 2 ) T AX 2 .j , j = 1 — >NX 2 (7) 

The above equations state mathematically the fun- 
damentally important concept that in a system op- 
timization the contributing disciplines should not 
optimize themselves for a traditional, discipline- 
specific objective such as the minimum aerodynamic 
drag or minimum structural weight. They should 
optimize themselves for a “synthetic” objective 
function that measures the influence of the BB r de- 
sign variables X r on the entire system objective 
function. 

Another way to look at it is to observe that, in long- 
hand 

<l> 2 = D(y,,„ x 2J ) t AX 2 ,i + D(y u , x 2 . 2 ) t AX 2i2 +... 

+ D(yi,j, x 2 j) T AX 2 j + .... j = 1 — >NX 2 (8) 

so it may be regarded as a composite objective function 
commonly used in multiobjective optimization. One 
may say, therefore, that in a coupled system the local 
disciplinary or subsystem optimizations should be 
multiobjective with a composite objective function. The 
composite objective should be a sum of the local design 
variables weighted by their influence on the single ob- 
jective of the whole system. It should be emphasized 
that this is true also in that particular BB r where O is 
being computed. In the aircraft example it is O = yj a in 
BB t according to eq. 1. However, the BBj optimization 
objective is not <(>, = y u . Instead, it is <J>i from an equa- 
tion analogous to eq.8. 

The local optimization problem may be stated formally 
for BB 2 


Given: 

X 2 , Z, and Y 2 ,, , Y 2>3 

(9) 

Find: 

AX 2 ; length NX 2 


Minimize: 

<l> 2 = D(y u ,X 2 ) T AX 2 


Satisfy: 

G 2 < 0, including side-constraints 


Incidentally, we adhere to the convention which calls 
for minimization of the objective function. If the appli- 
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cation requires that function be maximized, as it does 
in the example of aircraft range, we convert the objec- 
tive, e.g., O = - (range). 

The optimization problem for BBi , and BB 3 are anal- 
ogous. All three problems being independent of each 
other may be solved concurrently. This is an opportu- 
nity for concurrent engineering and parallel process- 
ing. 

By solving eq. 9 for all three BBs, we have improved- 
the system because, according to eq. 5 and 6, we have 
reduced $ by AO, while satisfying constraints in each 
BB. 

2.2. System-level optimization . 

So far we have improved the system by manipulating X 
in the presence of a constant Z. We can score further 
improvement by exploiting Zs as variables. To do so 
we need to know how Z influences O = yij. That is, we 
need D(y u ,Z). 

At this point, the BLISS algorithm forks into two al- 
ternatives, termed BLISS/A and BLISS/B. 

2.2.1. BLISS/A 

This version of BLISS computes the derivatives of Y 
with respect to Z by modified GSE, eq.(2.1/2) (equa- 
tions from other sections are cited in (), the section 
number given before the /). The GSE modification ac- 
counts for the fact that optimization of a BB turns its X 
into a function of Y and Z that enter that particular BB 
as parameters. The modification leads to a new gener- 
alization of GSE that takes the following form 



termed GSE/OS for GSE/Optimized Subsystems. The 
GSE/OS yields a vector D(Y,z k ) and D(X,z k ), and be- 
cause O is one of the elements of Y, d> = yi- ly we get the 
desired derivative D(d>,z k ). Derivation and details of 
the GSE/OS structure, including the definition of the 
matrix M, are in Section 2 of the Appendix. At this 
point it will suffice to say that the matrix of coefficients 
in GSE/OS is populated with d(Y r ,Y s ), d(Y r ,X r ), and 
d(X r ,Y s ). These terms and the RHS terms of d(Y,z k ) 
and d(X,z k ) are obtained from the following sources 

• d(Y r ,Y s ), d(Y r ,X r ), and d(Y,z k ) from BBSA 

• d(X r ,z k ), d(X r ,Y s ) from BBOSA 


The terms d(X,z k ) and d(X r ,Y s ) are the derivatives of 
optimum with respect to parameters that, in principle, 
may be obtained by differentiation of the Kuhn -Tucker 
conditions, e.g., an algorithm described in Sobieszc- 
zanski-Sobieski et al, 1982. That approach, however, 
requires second order derivatives of behavior, too 
costly in most large-scale applications. Therefore, an 
approximate algorithm adapted from Vanderplaats and 
Cai, 1986, is given in Section 3 of the Appendix. In 
that algorithm, parameters are perturbed by a small 
increment, one at a time, and the BB optimization is 
repeated by Linear Programming (LP) starting ffom 
the optimal point. Derivatives of optimal X and Y 
with respect to parameters are then computed by finite 
differences. 

2.2.2. BLISS/B 

This version of BLISS avoids calculation of d(X r ,z k ) 
and d(X r ,Y s ) altogether by using an algorithm that 
yields D(<I>,P), where P includes both Y and Z. The 
algorithm, described in literature (e.g., Barthelemy and 
Sobieszczanski-Sobieski, 1983) is based on the well- 
known notion that the Lagrange multipliers may be 
interpreted as the prices, stated in the units of O, for 
the constraint changes caused by incrementing p,. For a 
general case of the objective F=F(P) and G 0 = G 0 (P), the 
algorithm gives the following formula for D(F,P) C 

D(F,P) 0 =d(F,P) + L T d(G 0 ,P) 

To use the above in BLISS, consider that in P we have 
an independent Z but Y=Y(Z) so that the terms d() 
require chain -differentiation. Hence, the above general 
formula tranforms to 

D(y u ,Z) 0 T = (L t d(G 0 ,Z»] + (L T d(G„,Z)) 2 + 

(L T d(G ot Z)) 3 + [(L T d(G 0 ,Y))i+ (L T d(G 0 ,Y)) 2 
+ (L T d(G 0 ,Y)) 3 ](D(Y,Z)) + D(y u ,Z) T ( 1 ) 

where L is the vector of Lagrange multipliers and ( )i, 

( ) 2 , and ( ) 3 identify the BBs 1, 2, and 3. 

The terms in the above equation originate from the 
following sources: 

• d(G 0 ,Z) and d(G 0 ,Y) - BBSA performed on iso- 
lated BB r 

• L - obtained for BB r at the end of BBOPT 

• D(Y,Z) - from GSE in SSA 

• D(y u ,Z) - the column corresponding to y u in the 

above matrix D(Y,Z) 
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BLISS/B is substantially simpler in implementation 
than BLISS/A and it eliminates the computational cost 
of one LP per parameter Y and Z. Optimizers that yield 
L as a by-product of optimization are available for use 
in BBOPT, or L may be obtained as described in 
Haftka and Gurdal, 1992. 

2.2,3. Optimization in the Z-space. 

Once D(y K j,Z) have been computed from either 
eq.(2.2.1/l) or as D(y 1(1 ,Z) 0 from eq. (2.2.2/ 1), we can 
further improve the system objective by executing the 
following optimization, using any suitable optimizer 


sequence of these operations repeated with the new 
values of all quantities involved, including new values 
of all the derivatives because they would change if 
there were any nonlinearities in the system (as there 
usually are). 

In a large-scale application where execution of each 
BLISS cycle may require significant resources and 
time, the engineering team may wish to review the 
results before committing to the next cycle. That inter- 
vention may entail a problem reformulation, such as 
overriding the variable values, deleting and adding 
variables, constraints, and even BBs. 


Given: 

Z and d>o 

(1) 

Find: 

AZ 


Minimize: 

O = <J> 0 +D(y u ,Z) T AZ 


Satisfy: 

ZL < Z + AZ < ZU; AZL < AZ < AZU 


Where <I>o is inherited from the previous cycle SA for 
X and Z (initialized if it is the first cycle). It is rec- 
ommended to handle the Z constraints by means of a 
trust-region technique, e.g., Alexandrov 1996. In the 
above, term D(y ](i , Z) is a constrained derivative that 
protects G 0 = 0 in all BBs. Therefore, the optimization 
is unconstrained except of the side-constraints and 
move limits. 

However, some BBs may have constraints that depend 
on Z and Y more strongly than on X (in the extreme 
case some constraints may not be functions of X at all, 
only of Y and Z). Such constraints, denoted G yz , may 
be difficult (or impossible) to satisfy by manipulating 
only X in BBOPT. To satisfy them, one must add to the 
Z-space optimization in eq.l their extrapolated values 

G yz T =G yz . 0 T + (d(G yz ,Z) + 

d(G yz ,Y)D(Y,Z)) T AZ < 0 (2) 

where d(G yz ,Z), and d(G yz ,Y) are obtained from BBSA. 
In this instance, the Z-space optimization becomes a 
constrained one. 

3. Iterative Procedure 

The two operations, the local optimizations in the BBs 
and the system-level optimization, described in Sec. 2, 
result in a new system, altered because of the incre- 
ments on X and Z. This means that inputs to and out- 
puts from SA, BBA, BBSA, SSA, BBOPT, BBOSA 
(BLISS/ A), and SOPT all need to be updated, and the 



Figure 2: BLISS/B Flowchart 


Thus, the following procedure emerges, illustrated also 
by a flowchart in Figure 2 for BLISS/B with the 
BLISS/A operations, if different, noted in f ]. 

0. Initialize X & Z. 

1. SA to get Ys and Gs; this includes BBAs for all 
BBs. 

2. Examine TERMINATION CRITERIA, exercise 
judgment to override the results, modify the problem 
formulation, and CONTINUE or STOP. 

3. BBSA to obtain d(Y,X), d(Y r . s ,Y s ), d(G,Z), and 
d(G,Y), and SSA, eq. (2.1/2), to get D(Y,X) [and 
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D(Y,Z)]; Here is an opportunity for concurrent process- 
ing . 

4. BBOPT for all BBs, eq. (2.1/9) using 
formulatedindividually for each BB (eq. (2.1/6, 7)), 

get (j)opt and AXo pt ; obtain L for G 0 [skip L]. Here is an 
o pportunity for concurrent processing . 

5. Obtain D(0,Z) as in eq. (2. 2.2/1). [Execute BBOSA 
to obtain d(X,Z) and d(X,Y), and form and solve 
GSE/OS (Appendix, Section 3) to generate D(Y,Z)]. 
Here is an opportunity for concurrent processing . 

6. SOPT to get AZopt by eq. (2. 2.3/1 and 2) herein. 

7. Update all quantities, and repeat from 1 . 

X = X 0 +AX opt ;Z = Z + AZ opl 

Note: Termination is placed as #2 after SA to ensure 
that the full analysis results document the final system 
design, as opposed to having it documented only by the 
extrapolated quantities. Also, at this point the engi- 
neering team may decide whether to intervene by 
modifying the variable values, and adding or deleting 
the design variables and constraints. 

When started from a feasible design, the procedure will 
result in an improved system, while the local con- 
straints are kept satisfied within extrapolation accu- 
racy, even when terminated before convergence. 



Figure 3: Polynomial Representation of Wing Twist 

In case of an infeasible design start, the improvement 
will be in the sense of reductions in the constraint vio- 
lations, while the objective may exhibit an increase, at 
least initially. The procedure achieves the improvement 
by virtue of optimization alternating between the do- 


main of NB X-spaces (Step #4) and the single Z-space 
(Step #6). 

Caveat: because in BLISS/B the extrapolation of 3> in 
eq. (2.2.3/1) is based on the Lagrange multipliers in 
eq. (2.2.2/1), its accuracy depends on the BBOPT 
yielding a feasible solution, and on the active con- 
straints G 0 remaining active for updated Z. If some 
constraints leave the active set G 0 , or new constraints 
enter, a discontinuous change of the extrapolation error 
may result. For example, consider the wing aspect ratio 
AR as a Z-variable and suppose that for AR = 3 it is 
the stress due to the wing bending that is one of the 
active constraints in the structures BB. If optimization 
in the Z-space took the design to AR = 4, the next cy- 
cle may reveal that the stress constraint is satisfied but 
a flutter constraint becomes critical. Past experience 
(Sobieszczanski-Sobieski, 1983) shows that this dis- 
continuity is likely to slow, but not to prevent, the pro- 
cess convergence, and may be controlled by adjusting 
the move limits. 



Figure 4: Data Dependencies for Range Optimization 

4. Numerical Tests and Examples 

BLISS/ A was tested on a sample of test problems from 
Hock and Schittkowski, 1981, and on a design of an 
electronic package. BLISS/B was exercised on the lat- 
ter, and also on a very simplified aircraft configuration 
problem. Both versions of BLISS performed as in- 
tended in all of the tests. The sole purpose of these 
initial numerical experiments was to test and to dem- 
onstrate the BLISS procedure logic and data flow, 
therefore, the BBs were merely surrogates of the nu- 
merical processes that need to be used in real applica- 
tions. 
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4.1. Aircraft Optimization 

The aircraft test was an optimum cruise segment of a 
supersonic business jet based on the 1995-96 AIAA 
Student Competition. This problem was selected be- 
cause of its available data base and the availability of 
the black boxes written in Visual Basic in form of Ex- 
cel spreadsheets. The supersonic business jet was mod- 
eled as a coupled system of structures (BBj), aerody- 
namics (BB 2 ), propulsion (BB 3 ), and aircraft range 
(BB 4 ). All the disciplines were represented by modules 
comprising an analysis level typical for an early con- 
ceptual design stage. 


var \ cycle* 

1 

2 

3 

4 

5 

Range (SSA) 

535.79 

1581 .67 

3425.35 

3961.41 

3963.98 

Extpl. Error 

-535.79 

-536.67 

-431.63 

-56.26 



■BQ 

-0.16 

-3.26 

-0.86 




0.00 

0.00 

0.00 

0.00 

■ in jinw 



HHZXH 


0.00 

X Extpl. 



■Ezra 

-0.86 

0.00 

Z Extpl. 

449.19 

1301.30 



0.00 

Range (Extpl.) 

■utram 

2993.72 

3905.15 

■sHdMdd 


X 

0.25 

0.14951 

0.17476 

0.25775 


X 

1 

0.75 

0.75 

0.75 

0.75 

c, 

1 

0.75 

0.75 

0.75 

0.75 

T 

0.5 

0.1676 

0.20703 

0.15624 

0.15624 

t/c 

0.05 

0.06 

0.06 

0.06 

0.06 

h(ft) 

45000 

54000 

60000 

60000 

60000 

M 

1.6 

1.4 

1.4 

1 .4 

1 .4 

AR 

5.5 

4.4 

3.3 

2.5 

2.5 

A(°) 

55 

66 

70 

70 

70 


1000 

1200 

1400 

1500 

1500 


*One cycle is one pass through the BLISS procedure 


Table 1: A/C Results for 20% Move Limit 

The aircraft optimization was a maximization of ther- 
ange computed through the Breguet range equation. 
For testing purposes, additional design and state vari- 
ables were introduced in BBs 1 through 3, and func- 
tional relationships not present in the original BBs 
were supplied to reflect what is commonly known 
about the typical functions involved in design. For ex- 
ample, stress is expected to fall as a reciprocal of the 
increase of the skin thickness in a wing box. Such re- 
lationships were represented by polynomial functions. 
One plot of such a function is shown in Figure 3, por- 
traying the wing twist as a function of the wing box 
cross-sectional dimensions scale factor and the wing 
lift. 

Section 4 of the Appendix defines the BBs in this ex 
ample by their input and output variables, and by the 
functions that link output to input. Table A1 also iden- 
tifies local constraints and side constraints. Note that 
BB 2 contains a constraint that does not depend on its X 
or Y input, thus the Z-space optimization is a con- 
strained one, per eq. (2.2.3/1 & 2 ). Side constraints on 


Z were judiciously selected to guard against conditions 
not accounted for in the BBAs. For example, the lower 
bound of 2.5 on aspect ratio stemmed from the sub- 
sonic performance considerations. 


num \ den 

X 

X 

C, 

T 

tfc 

w, 

0.01 146 

1.71536 

0.01981 

-0.15744 

0.12714 

W F 

0 

0 

0 

0 

0.72626 

e 

-0.03342 

0.19971 

3.31E-15 

-1.73E-14 

-2.10E-14 

L 

0.01 146 

1.71536 

001981 

-0.15744 

0.12714 

O 

-4.19E-05 

0.00581 

0.12457 

-0.00049 

0.68106 

LTD 

0.0115 

1.7095 

-0.1046 

-0.15694 

-0.54935 

SFC 

1.98E-20 

-5.07E-18 

-2. 70 E- 17 

0.08544 

0 

W E 

-4.40E-05 

0.0061 

0.13003 

-1.03986 

0.71531 

ESF 

-4.19E-05 

0.00501 

0.12457 

-0.99059 

0.68108 

R 

-0.00077 

-0.12692 

-0.12581 

-0.07299 

0.10115 

num \ den 

h 

M 

AR 


S,* 

W T 

-0.33931 

0.31958 

0.08208 

0.2537 

■EE5E3 

W F 




IEHK] 

■HE3D 

e 

BB213U 



k zjzaa 

K3E2Z11 

L 


0.31958 

■3KI 

MWflPESTlI 


D 






UD 






SFC 



Exilian 

BE3EB 


w e 


MREM 

■ft&fsEra 

BBEEBSI 

1.04857 

ESF 

■ 

BEE5E53I 


MgjEEESII 

0.99838 

R 

2.076161 

-1.04784 

-0.39618 

0.82904 

0.15535 


Table 2: Normalized Y Derivatives w.r.t. X and Z 


The BBs are coupled by the output-to-input data 
transfers (design structure matrix) depicted in Figure 4. 
Note that BB 4 is an analysis-only module and does not 
feedback any data to other BBs. 



Figure 5: Range and Extrapolation Error Histogram 

This test was conducted entirely using MATLAB 5 and 
its Optimization Toolbox. The entire MATLAB code 
listing for the aircraft range model may be found in 
Section 5 of the Appendix. To establish a benchmark, 
the system was first optimized using an all-in-one ap- 
proach in which the MATLAB optimizer was coupled 
directly to SA and saw no distinction between the X 
and Z variables. Next, the test case was executed using 
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BLISS/B, starting at different infeasible initial points 
chosen by varying the six design variables that are not 
arguments in the polynomial functions. The choice of 
initial values for variables that are arguments of the 
polynomial functions was limited due to the nature of 
the polynomial formulation. This limitation is not a 
characteristic of the BLISS method itself, as the poly- 
nomial functions would not be required in a large scale 
optimization problem. With the move limits ranging 
from 10 to 70 %, the procedure convergence was sat- 
isfactory through the move limits of 60% for all initial 
points tested. However, in nearly all cases, no addi- 
tional improvement in convergence rate was recorded 
for move limits greater than 20%. For instance, the 
objective function was advanced to within 1% of the 
benchmark in 5 passes for move limits 20 and 30%. 
Onset of an erratic behavior was observed with move 
limits increased past 60%, the procedure converged or 
diverged dependent on the starting point. 



Figure 6: Range Sensitivities (1 st cycle) 

Table 1 displays a sample of typical results for the 
move limits value of 20%. It shows that the initial de- 
sign range was extremely poor, only 536 nm. BLISS/B 
improvements advanced the range to 3964 nm. The 
range converged monotonically, although in some 
cases small amplitude oscillations were observed. 
Comparison of the extrapolated and actual values of 
the objective and constraints shows reasonable accu- 
racy and conservatism of the extrapolations. The opti- 
mal values of the design variables reflect numerous 
tradeoffs typical for aircraft design. For instance, opti- 
mal t/c resulted, in part, from a trade-off between the 
wave drag and structural weight. Table 2 shows nor- 
malized (logarithmic) derivatives of all Ys , including 
the range, w.r.t. all the X and Z variables, sampled in 
Cycle 1 to illustrate sensitivity of the system solution to 
design variables. 


Figure 5 illustrates the range histogram, and depicts 
the extrapolation error as being effectively controlled 
by the move limits. Range sensitivities to X and Z 
variables are shown in Figure 6. As expected, altitude 
and Mach number have the largest effect on the objec- 
tive function, while taper ratio has the smallest. 



Figure 7: BB and System Contributions to Range 

Figure 7 shows the individual BB and system contri- 
butions to the range objective in each cycle. Here it is 
observed that, in this particular case, the contribution 
of system level variables is significantly larger than 
that of the local variables in the extrapolation of range. 

This test case was also implemented in a software 
package for system analysis and optimization called 
iSIGHT (iSIGHT, 1998). The iSIGHT and MATLAB 
results cross-check was completely satisfactory. 

4.2 Electronic Package optimization 

The electronic packaging was introduced as an MDO 
problem in Renaud, 1993. Its electrical and thermal 
subsystems are coupled because component resistance 
is influenced by operating temperatures and the tem- 
peratures depend on resistance. 

The objective of the problem is to maximize the watt 
density for the electronic package subject to con- 
straints. The constraints require the operation tem- 
peratures for the resistors to be below a threshold tem- 
perature and the current through the two resistors to be 
equal. The system diagram in Figure 8 shows the data 
dependencies for two BBs, representing electrical re- 
sistance analysis and thermal analysis. As Figure 8 
indicates, there are no “natural” Z’s in this case. 
Therefore, Z’s were created as targets imposed on each 
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of the Y’s and the BBOPTs were required to match 
the Y values to those Z targets (similar as it is done in 
the Collaborative Optimization method). Details of the 
electronic packaging problem may be found in Padula, 
1996. 


(Xi,Xi,Xi,X *) 


(X 5 ,X6,X7,X*) 


Re ristanre 
(Yz,Y 3 ) 


Temperature 

(Yn,Yu) 





/ 1 





Thermal 




Electrical 




. 

1 1/ ' 



) 


Figure 8: Electronic Packaging Data Dependencies 


This test case was implemented in iSIGHT using 
BLISS/A and B. A benchmark result was obtained by 
executing an all-in-one optimization from various 
starting points (“A-in-O” column). BLISS/A and B 
were started from the same points. Table 4 displays the 
benchmark and the BLISS/A and B results as showing 
a good agreement. Table 4 also indicates a comparison 
of the computational labor (the “Work” column) meas- 
ured by the number of BB evaluations necessary to 
converge the fixed-point iterations in BBAs and in SA, 
all repeated as needed to compute derivatives by finite- 
differences in a gradient-guided optimization. As Table 
4 shows, the BLISS/B computational labor was sub- 
stantially lower than the benchmark in all cases. 


Method 

Case 

Initial Design 
Objective 

knit. Dee. Max 
Conetr. Viol. 

Final Design 
Obj active 

Fin. Des. Max 
Conetr. Viol. 

Wort 

A-in-O 

1 

7.79440E+01 

2.16630E-08 

6.39720E+05 

1 -22E-03 

496 


2 

6.83630E+03 

-2.89560E-01 

6.39720E+05 

1.22E-03 

264 


3 

1.511 10E+03 

-4.292406-02 

6.36540E+05 

1.45E-03 

264 


4 

1.46070E+01 

-1.02490E-03 

6.36940E+05 

1 .42E-03 

175 

BLBS/A 

1 

7.7944OE+01 

2.16630E-08 

6.39700E+05 

1 .20E-03 

436 


2 

6.83630E+03 

-2.89560E-01 

6.39050E+05 

1.18E-03 

508 


3 

1.511 10E+03 

-4.29240E-02 

6.390 50E+05 

-4.896-04 

174 


4 

1 .46070E+01 

-1.024906-03 

6.39290E+05 

3.70E-04 

313 

BLISS/B 

1 

7.79440E+01 

2.16630E-08 

6.39720E+05 

1.22E-03 

365 


2 

6. 83630 E +03 

-2.89560E-01 

6.39720E+05 

1 .22E-03 

207 


3 

1.511 10E+03 

-4.29240E-02 

6.397206+05 

1 .22E-03 

114 


4 

1 .46070E+01 

-1.02490E-03 

6.397206+05 

1 .22E-03 

105 


Table 4: Electronic Packaging Data 

5. BLISS Status, Assessment, and Concluding 
Remarks 


A method for engineering system optimization was 
developed to decompose the problem into a set of local 
optimizations (large number of detailed design vari- 
ables) and a system-level optimization (small number 


of global design variables). Optimum sensitivity data 
link the subsystem and system level optimizations. 
There are two variants of the method, BLISS/A and 
BLISS/B, that differ by the details of that linkage. In 
the paper, the method algorithm was laid out in detail 
for a system of three subdomains (modules). Its gener- 
alization to NB subdomains is straightforward. The 
same algorithm may be used to decompose any of the 
local optimizations, hence optimization may be con- 
ducted at more than two levels. 

MATLAB and iSIGHT programming languages were 
used to implement and test the method prototype on a 
simplified, conceptual level supersonic business jet 
design, and a detailed design of an electronic device. 
Dimensionality and complexity of the preliminary test 
cases was intentionally kept very low for an expedi- 
tious assessment of the method potential before more 
resources are invested in further development. Favor- 
able agreement with the benchmark results and a sat- 
isfactory convergence observed in the above tests pro- 
vided motivation for such development and future 
testing in larger applications. 

Assessment of BLISS at the above development status 
is as follows. BLISS relies on linearization of a gener- 
ally nonlinear optimization, therefore its effectiveness 
depends on the degree of nonlinearity. As any gradi- 
ent-guided method, it guarantees a cycle-to-cycle im- 
provement, but if the problem is non-convex, its con- 
vergence to the global optimum depends on the starting 
point and may strongly depend on the move limits. In 
this regard, BLISS’s strong points are in the procedure 
being open to human intervention between the cycles 
and in the autonomy of the subdomain optimizations in 
local variables. These optimizations may be conducted 
by any means deemed to be most suitable by discipli- 
nary experts, hence non-convexity, and strong nonline- 
arities in terms of the local variables often encountered 
in subdomains, e.g., the local buckling in thin-walled 
structures, are isolated and prevented from slowing 
down the system-level optimization convergence. On 
the other hand, the optimization robustness may be 
adversely affected by the local constraints leaving and 
entering the active constraint set. Effect of the above 
on BLISS/A is much less than on BLISS/B. This is 
probably the only reason to continue the development 
of BLISS/A alongside with BLISS/B, even though 
BLISS/B has a distinct advantage in simplicity and a 
much lower computational cost. Once there is more 
information on the relative merits and demerits of both 
variants, the better variant may be selected. 

The demand BLISS puts on the computer storage is the 
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same the subdomains would require for their own, 
stand-alone optimizations, with exception of the gen- 
eration and solution of the Global Sensitivity Equa- 
tions. If there is a pair of BBs that exchange large 
number of the y riS>i - quantities, dimensionality of the 
corresponding matrices that store the derivatives, and 
computational cost of these derivatives needed to form 
GSE, may become prohibitive. Some relief may be pro- 
vided here by application of condensation techniques 
and by deleting from GSE those derivative matrices 
that are known to have negligible effect on the system 
behavior. 

The principal advantage of BLISS appears to lie in its 
separating overall system design considerations from 
the considerations of the detail. This makes the re- 
sulting mapping of its algorithm fit well on diverse, 
and potentially dispersed, human organizations. This 
advantage remains to be demonstrated in further de- 
velopment toward large-scale, complex applications. 
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Appendix 

This Appendix provides details of the Global Sensitiv- 
ity Equations (GSE) applied to a system which opti- 
mizes BBs, the details of a technique for the BB Opti- 
mum Sensitivity Analysis, and the details of the air- 
craft range optimization model. 


1. Global Sensitivity Equations 

Derivatives of Y w.r.t. X, and Z, are obtained rigor- 
ously from the Implicit Function Theorem in Sobieszc- 
zanski-Sobieski, 1990. The condensed derivation is 
provided below. It begins by recognizing that 


(Al) 

Y,.2 

= Y,.2(Z, 

X 2 , 

,Y 2 , 

> y 2 , 3 ) 

(A2) 

Yu 

= Yu(Z, 

X 3 , 

Y 3 .! 

. Y3.2) 

(A3) 

y 2 .. 

N 
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X,, 

Y,,2 

. Yu) 

(A4) 

Y 2 ,3 

= Y 2 .3(Z, 

X3, 

Y3.1 

. Y3.2) 

(A5) 

y 3 ., 

= Y 3 ,,(Z, 
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Y,,2 

, Yu) 

(A6) 

Y3.2 

<N 

tn 

> 

II 

X 2 . 

y 2iI 

. Y 2i3 ) 


where the independent variables are X and Z. 

Observe that eq. A1-A6 are coupled by Y it e.g., Y 3t3 
depends on Y u in eq. A5, and Y I>3 depends on Y 3>1 in 
eq. A2. Consider for an example, the chain- 
differentiation w.r.t. xij applied to eq. A3. It yields 

A7) D(Y 2 , 1 ,x 1 , j ) = d(Y 2J ,x 1J ) + 

d(Y 2i i,Yi)D(Y 1 ,x ltj ) + 

d(Y 2 ,i,Y 3 ) D(Y 3 ,X]j) 

Repeating the above for the remaining equations, 
treating Y 2 ,i as a subset of Yj, and collecting the terms 
leads to eq. (2.1/2 and 3). 

The derivatives of Y w.r.t. z k are obtained by simply 
replacing x rJ with z k in eq. (2.1/2) to obtain 

A8) [A]{D(Y,z k )} = {d(Y,z k )} 

2. GSE/Optimized Subsystems 

In the preceding section both X and Z are independent 
variables. By virtue of BBOPT conducted for constant 
Z and Y inputs, X becomes dependent on Z so that 
derivatives of X w.r.t. exist in addition to derivatives of 
Y w.r.t. Z. For example, optimal X 2 depends on Z, 
Y 2 j i, and Y 23 , that are parameters in the optimization 
of BB 2 . Hence, to compute the derivatives of Y and X 
w.r.t. Z, we begin by rewriting the functional relation- 
ships in eq. A1-A6, adding the new dependencies in all 
three BBs in the system, 

(A9) Y 1 , 2 = Y 1 , 2 (Z.X 2 „Y 2 , 1 ,Y 2 , 3 ) 
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(A 10) 

Y i ,3 = Y 1 , 3 (Z. x 3 , y 3 ,„ y 3 . 2 ) 

derivative expressions is written it forms a set of si- 
multaneous, algebraic equations in which the total de- 

(All) 

Yj., - Y 2 .,(Z, X u Y ,.2* Y,. 3 ) 

rivatives such as D(Y 2 ,z k ) and D(X 2 ,z k ) appear as un- 
knowns. This is a new generalization of GSE, 

(A 12) 

Y 2 , 3 = Y 2 . 3 (Z, x 3 , y 3>1 ,y 3 . 2 ) 

termed GSE/OS for GSE/Optimized Subsystems. 

For the case of three-BB system, these equations may 

(A13) 

Y 3 .,=Y 3 .,(Z, X,. Y,. 2l Y u ) 

be presented in a matrix format like this 

(A 14) 

Y 3 . 2 = Y 3 , 2 (Z, X 2 , Y 2 .,.Y 2 . 3 ) 

(A20) [M yy ] { D(Y,z k ) } + [M yx ]{D(X,z k )} = d(Y,z k ) 
[M xy ]{D(Y,z k )} + [M xx ]{D(X,z k )} = d(X,z k ) 

(A15) 

X, = X)(Z, Y]. 2 , Y, i3 ) 

The internal structure of the M-matrices in the above is 

(A 16) 

X 2 =X 2 (Z, Yjj.Yjj) 

for [M yy ] : 

(A 17) 

X 3 = X 3 (Z, Y 3 .,. Y 3i2 ) 



The same Implicit Function Theorem that is the basis 
of the GSE derivation may be applied to the above 
equations to obtain D(Y,Z). For example, by applying 
chain-differentiation to Y 2 J treated as a subset of Y 2 , 
we obtain 

(A18) D(Y 2 ,z k ) = d(Y 2 ,z k ) + d(Y 2 ,X 2 )D(X 2 ,z k ) + 
d(Y 2 ,Y 1 )D(Y,.z k ) + 

d(Y 2 ,Y 3 )D(Y 3 ,z k ) 


and for X 2 , again as one example: 

(A 19) D(X 2 , z k ) = d(X 2 ,z k ) + d(X 2 ,Y l )D(Y lt z k ) + 

d(X 2 ,Y 3 )D(Y 3 ,z k ) 


In the above, the D-terms are the total derivatives we 
seek, while the d-terms are the partial derivatives of 
two, different kinds. The derivatives of Y r w.r.t. Y s and 
Y r w.r.t. X r are obtained from BBSA r using any sensi- 
tivity analysis algorithm appropriate for the particular 
BB r (including the option of finite differencing). The 
derivatives of X r w.r.t. z k and X r w.r.t. Y s are produced 
by an analysis of optimum for sensitivity to parameters, 
BBOSA r , explained in later in this Appendix. 

As a mathematical digression, one should mention at 
this point that the derivatives termed partial in the 
above would be called total in both BBS A and BBOSA. 
This is not a contradiction. It is so because the partial 
and total derivatives are hierarchically related in a 
multilevel system of parents and children. What is a 
total derivative in a child is partial at the parent level. 
In the application herein, the system of coupled three 
BBs is a parent, each BB is a child. 

The chain-derivative expressions for Yj, Y 3 , Xi and X 3 
look similar to eq. A18 and A 19, differences are only 
in the subscripts. When the entire set of six chain - 


for [M yx ] : 


for [M xy ] : 


and for [M xx ] : 



Again, in the above, all Y rs are folded into Y r for 
compactness, and the terms are falling into the previ- 
ously introduced categories as follows: 

• M yy , M yx , and d(Y,z k ) from BBS A 

• M xy and d(X,z k ) — from BBOSA 

As in GSE, one may obtain D(Y 2 ,z k ) and D(X 2 ,z k ) for 
all z k , k = 1— >NZ by means of one of the efficient 
techniques for linear equations with many right-hand- 
sides. 


3. Black Box Optimum Sensitivity Analysis (BBOSA) 
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Analysis of optimum for sensitivity to parameters (also 
called the post-optimum analysis) is preceded by solv- 
ing a BB optimization problem 


(A21) Given: 
Find: 

Minimize: 

Satisfy: 


P 

X 

(X,P) 

G(X,P) < 0, including side- 
constraints and move limits 


z k , and Y r>s . because these quantities are kept constant 
in BBOPT r . 

After <) ) min , and Xo pt are found, one may seek sensitivity 
of these quantities to the change of P in form of the 
derivatives D^^P) and D(X opl ,P). 

Vanderplaats and Cai, 1986, review techniques, rigor- 
ous and approximate, available for calculating 
D(Xo pt ,P). The technique adapted for the BLISS/A 
purposes comprises the following steps executed for 
BB r : 


where P are parameters kept constant while an opti- 1. Choose parameter P k , an element of Z or Y, and in- 
mizer manipulates X. crement it by a AP 


In the BLISS application, the parameters P in BB r are 2. Use derivatives from SSA to extrapolate F and G 0 


Table Al: BB Definitions 


BB 

EJP99 

Internal 


Structures 

AR, A, 

t/ Q 
/C ’^REF* 

W P0 ,W 0 , 
W E , L, 
N z ,X,x 

t c 

C ° REF 1 1 + Tk 

t c ; b/2-JS REF AR /2; R= ; 0 = 

^ befAR 3(1 + X) 

pf(x,^/,R,L); Fol = pf(x); W w = (0.005 1(W T N z ) ft557 S REF °‘ w 

AR “(^/) “"(I + X )° '(0.1 875S ref ) 01 / cos ( A ))Fo 1 ; W FW = 

(5S,„/18)(% 1)42.5; W F = W FW + W FO ; W T = W„ + W w + W F 

+ W E ; o 1 -» a 5 = pf(^/,L, x, Jj^.R); 

W T ,W F , 

e 

Constraints! al -> a5 < 1.09; 0.96 < © < 1.04 

Aero- 

dynamics 

M,h, A, 

ar ,>{, 

S W 

J REF » ” T* 

0 , ESF, 

r 

^ Dmin. M<l» 

c, 

if h < 36089ft, V = Ml 1 16.39^1 - (6.875e - 06)h , p = (2.377e - 03)( 1 
-(6.875e - 06h)) 42561 ; V = M968.1, p = (2.377e - 03)e (h 360891 20 8 06 7 , 
W T 

if h > 36089 ft; C L = } ; Fol = pf(ESF,C,); C n = 

0.5pV 2 S REF ' D ““ 

C D„i.. M <,Fol+3.05()/) S5 cos(A) 5j ; k = 1 / (it0.8AR); Fo2=pf(0); 
C d = ( C d„ 10 + kC L ’)Fo2; L = W t ; D = C D 0.5pV 2 S REF ; 
dp/dx = pf(//) 

r 

a 

dp/dx ^ 1 .04, evaluated at system level 

Propulsion 

M, h, D, 

W be.T 

T = T *161 68.6; Temp = pf(M,h, T); ESF = (^/)/T; 

SFC = 1.1324 + 1.5344M - (3.295 6e - 05)h - (1.637 9e - 04 )T 
-031623M 2 + (821 38e- 06)Mh - ( 10.496e - 05)TM - (8.574e - 1 l)h 2 
+ (3.8042e- 09)Th + (1 ,0600e - 08)T 2 ; W E = 3W be ESF 105 ; 

Tim = 1 1484+ 10856M - 0.50802h + 3200.2M 2 - 0.29326Mh 
+(6.8572e - 06)h 2 

SFC, W F , 
ESF 

Constraints) 0.5 S ESF S 1.5; TsTm; Temp < 1.02 

Range 

M, h, L . 
D 

w T ,w F , 

SFC 

6 = 1 - 6.875e - 06* h, if h < 36089ft; 0 = 0.7519 

if h> 36089ft; R= M(LD)661 ^ ta f W T ] 

SFC W T - W F 

V * F J 

R 

Constants 


loom. W„- ZiOOMte M,.6g w„.«60lb; C„ -0.01 


Side Con- 
straints 

0.1 5 X < 0.4; 0.75 < x £ 1 .25; 0.75 s C f S 1.25; O.UTsl.0; 0.01 5 y < 0.09; 
30000 < h < 60000; 1.4<M<1.8; 2.5 < AR ^ 8.5; 40 < A < 70; 500 < S RFF < 1 500 
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linearly and by Linear Programming solve 

Find X 

Minimize F 
Satisfy G 0 < 0 

XL <= X <= XU; where XL and XU incorporate 
the side constraints and the move limits; 

to obtain X opt . 

3. Approximate D(X,P) = AX opt /AP 

4. Repeat from #1 for all elements of Z and Y input 
into BB r . 

Repeated for all BBs, the above procedure yields a set 
of D(X,Z) and D(X,Y) to be entered as d(X,Z) and 
d(X,Y) into GSE/OS, eq. A20. Solution of eq. A20 
provides D(X,Z) and D(Y,Z). The latter is substituted 
into eq. (2. 2.3/2), and D(0,Z), extracted from D(Y,Z), 
goes into eq. (2.2.3/1). 

4. A/C Range Optimization Model . 

Table A1 shows the equations used in each of the BBs 
for the aircraft model. Polynomial functions are repre- 
sented by ‘pf()’ with independent variables in the pa- 
rentheses. Each polynomial function is of the form: 

(A22) PF = A„ + A,*S T + (l/2)*S*A ij *S T 

Where S is the vector of independent variables, and A*, 
Ai, and Ay are coefficient terms. 

In calculating the polynomial functions using eq. A22, 
terms in the S vectors are in the same order as they 
appear in pf() in Table Al. The off diagonal terms of 
Ay are random numbers between 0 and 1. For this 
model, they are 


-- 

0.3970 

0.8152 

0.9230 

0.1108 

0.4252 

-- 

0.6357 

0.7435 

0.1138 

0.0329 

0.8856 

__ 

0.3657 

0.0019 

0.0878 

0.7248 

0.1978 

-- 

0.0169 

0.8955 

0.4568 

0.8075 

0.9239 

-- 


The remaining coefficient are: 

• 0 — > Ao = [1.0]; A = [0.3 -0.3 -0.3 -0.2]; 

Ay = [0.4 -0.4 -0.4 0]; 

• Fol — > A, = [1.0]; A, = [6.25]; Ay = [0]; 

• ol — > Ao = [1.0]; A = [-0.75 0.5 -0.75 0.5 






0.5]; Ay =[-2.5 0-2.5 0 0]; 
o2 — > Ao = [1.0]; A ; = [-0.5 0.333 -0.5 0.333 

0.333]; Ay = [-1.111 0-1.111 0 0]; 
o3— > Ao = [1.0]; A = [-0.375 0.25 -0.375 

0.25 0.25]; Ay = [-0.625 0 ^0.625 0 0]; 
g 4 — > Ao = [1.0]; A = [-0.3 0.2 -0.3 0.2 0.2]; 

A u = [-0.4 0-0.4 0 0]; 

g 5 — > Ao = [1.0]; Ai = [-0.25 0.1667 -0.25 

0.1667 0.1667]; Ay = [-0.2778 0 
-0.2778 0 0]; 

Fo2 — > Ao = [1.0]; A, = [0.2 0.2]; Ay = [0 0]; 
Fo3 — > Ao = [1.0]; Ai = [0]; Ay = [0.04]; 
dp/dx — > Ao = [1.0]; A, = [0.2]; Ay = [0]; 

Temp — > Ao = [1.0]; A, = [0.3 -0.3 0.3]; 

Ay = [0.4 -0.4 0.4]; 


Equations for SFC and the upper constraint bound on 
throttle setting in the Propulsion BB are polynomials 
representing surfaces fit to engine deck data (AIAA/ 
UTC/Pratt & Whitney, 1995/96). 


5. A/C Range MATLAB Code . 

Included in the following pages is the MATLAB code 
for the aircraft range optimization model. The con- 
strained optimization routine used in BBIOPT, 
BB20PT, BB30PT, and SYSOPT may be found in 
MATLAB’s Optimization Toolbox and is based on a 
Sequential Quadratic Programming method. The finite 
differencing subfunctions in FIN_DIFF are simple one- 
step forward finite difference codes that use a 1 percent 
step increment. 
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% them. Local optimizations are performed on each BB (BBOPT) as well as a 

% system level optimization (SOFT) using a gradient guided path based on the 
ogram Listing % Lagrange multipliers (OSAAA). Finally, all optimized changes to design 

% variables are used to update the model for an improved range. 

Page Number % 
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GRADphi4_Z = [Dphil_Z + Dphi2_Z + Dphi3_Z + Dphil_Y*DYEl_Z + 
DY_tc=A\dY_tc; Dphi2_Y*DYE2_Z + Dphi3_Y*DYE3_Z + DY4_Z] 
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% Subfimctions 

Subfunction BB2WRAPPER % BB_dragpolar -Calculates aerodynamic values 
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